#include "effective_ensemble.hpp"
#include "transition_map.hpp"
#include "channels.hpp"
#include "algorithm.hpp"
#include "input.hpp"
#include <boost/numeric/ublas/io.hpp>
#include <boost/bimap.hpp>
#include <fstream>

using namespace boost::numeric::ublas;
int main(int argc,char* argv[]){

/////////////////////////////////////////////////
  // Warnings if program is not called correctly //
  /////////////////////////////////////////////////
  if(argc<2){
    std::cerr << "Input file must be specified"<<std::endl;
    exit(1);
  }
  if(argc>2){
    std::cerr << "Too many arguments!"<<std::endl;
    exit(1);
  }
  /////////////////////
  // Read input file //
  /////////////////////
  std::string input_fn=std::string(argv[1]);
  input in(input_fn);
  ////////////////
  // Initialize //
  ////////////////
  int N=in.get_sample_size();
  double T=in.get_time_length();
  double dt=in.get_timestep();
  int dim_of_H=in.get_dim_of_H();
  vector<double> P0=in.get_initial_distribution();
  matrix<std::complex<double> > states0=in.get_initial_states();
  matrix<double> params=in.get_parameters();
  matrix<double>  J=in.get_operators();
  int number_of_channels=in.get_number_of_channels();
  
  std::string ens_out=in.get_ens_out();
  std::string P_out=in.get_P_out();
  std::string rho_out=in.get_rho_vec_out();
  std::string sigma_out=in.get_sigma_out();
  
  effective_ensemble eff(states0,P0);
  
  
  transition_map TLam[2];
  TLam[0].insert(position(0,1));
  TLam[1].insert(position(0,2));
  
  dtcl4 rates(params);
  
  channels chan(dim_of_H,number_of_channels,rates);
  chan.set_op(J);
 
  output out(3,2,N,ceil(T/dt),ens_out,P_out,rho_out,sigma_out);
  
  algorithm alg(dt,T,N,out);

  /////////////////
  // Calculation //
  /////////////////
  
  alg.apply(eff,chan,TLam);
  
  ////////////
  // Output //
  ////////////
  alg.print_to_files();
  
  return 0;
}


